Introduction to mrgsolve 2

Shen Cheng

2024-10-11

Outline

  • Model parameters
  • Model specifications
  • Model validations

An example model


[ prob ]
  ECP8506 mrgsolve model for illustration
  This model requires mrgsolve >= 1.0.3

[ pkmodel ] cmt = "GUT,CENT,PERIPH", depot = TRUE

[ param ]  
  WT   = 70
  EGFR = 90
  ALB  = 4.5
  AGE  = 35

[ theta ] @annotated
  0.443   : log(TVKA)
  4.12    : log(TVV2)
  1.17    : log(TVCL)
  4.21    : log(TVV3)
  1.28    : log(TVQ)
  0.485   : log(CLEGFR)
  -0.0377 : log(CLAGE)
  0.419   : log(CLALB)

[ omega ] @block 
  0.219
  0.0668 0.0824
  0.121  0.0703 0.114

[ sigma ] @block
  0.0399

[ pk ] 
  double V2WT   = log(WT/70.0);
  double CLWT   = log(WT/70.0)*0.75;
  double CLEGFR = log(EGFR/90.0)*THETA(6);
  double CLAGE  = log(AGE/35.0)*THETA(7);
  double V3WT   = log(WT/70.0);
  double QWT    = log(WT/70.0)*0.75;
  double CLALB  = log(ALB/4.5)*THETA(8);

  double KA  = exp(THETA(1) + ETA(1));
  double V2  = exp(THETA(2) + V2WT + ETA(2));
  double CL  = exp(THETA(3) + CLWT + CLEGFR + CLAGE + CLALB + ETA(3));
  double V3  = exp(THETA(4) + V3WT);
  double Q   = exp(THETA(5) + QWT);

  double S2 = V2/1000.0; //; dose in mcg, conc in mcg/mL

[ error ] 
  double IPRED = CENT/S2;
  double Y = IPRED * (1+EPS(1));

[ capture ] CL V2 IPRED Y

Load the model

suppressPackageStartupMessages(library(here)) 
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(mrgsolve))

mod <- mread(model = "106.cpp", 
             project = here("wk6/mrg-model"))
Building 106_cpp ... done.

Check the model

mod
param(mod)
omat(mod)
smat(mod)
see(mod)

Model parameters

mrgsolve maintains a parameter list (name / value pairs):

  • This list is used by default (under [ param ] or [ theta ]) if nothing else is done.
  • The parameter values in this list can be updated.

[ prob ]
  ECP8506 mrgsolve model for illustration
  This model requires mrgsolve >= 1.0.3

[ pkmodel ] cmt = "GUT,CENT,PERIPH", depot = TRUE

[ param ]  
  WT   = 70
  EGFR = 90
  ALB  = 4.5
  AGE  = 35

[ theta ] @annotated
  0.443   : log(TVKA)
  4.12    : log(TVV2)
  1.17    : log(TVCL)
  4.21    : log(TVV3)
  1.28    : log(TVQ)
  0.485   : log(CLEGFR)
  -0.0377 : log(CLAGE)
  0.419   : log(CLALB)

[ omega ] @block 
  0.219
  0.0668 0.0824
  0.121  0.0703 0.114

[ sigma ] @block
  0.0399

[ pk ] 
  double V2WT   = log(WT/70.0);
  double CLWT   = log(WT/70.0)*0.75;
  double CLEGFR = log(EGFR/90.0)*THETA(6);
  double CLAGE  = log(AGE/35.0)*THETA(7);
  double V3WT   = log(WT/70.0);
  double QWT    = log(WT/70.0)*0.75;
  double CLALB  = log(ALB/4.5)*THETA(8);

  double KA  = exp(THETA(1) + ETA(1));
  double V2  = exp(THETA(2) + V2WT + ETA(2));
  double CL  = exp(THETA(3) + CLWT + CLEGFR + CLAGE + CLALB + ETA(3));
  double V3  = exp(THETA(4) + V3WT);
  double Q   = exp(THETA(5) + QWT);

  double S2 = V2/1000.0; //; dose in mcg, conc in mcg/mL

[ error ] 
  double IPRED = CENT/S2;
  double Y = IPRED * (1+EPS(1));

[ capture ] CL V2 IPRED Y WT EGFR ALB AGE

Specify model parameters

[ param ]  
  WT   = 70
  EGFR = 90
  ALB  = 4.5
  AGE  = 35

[ theta ] @annotated
  0.443   : log(TVKA)
  4.12    : log(TVV2)
  1.17    : log(TVCL)
  4.21    : log(TVV3)
  1.28    : log(TVQ)
  0.485   : log(CLEGFR)
  -0.0377 : log(CLAGE)
  0.419   : log(CLALB)
[ param ]  
  WT   = 70
  EGFR = 90
  ALB  = 4.5
  AGE  = 35

[ theta ] 
  0.443   
  4.12    
  1.17    
  4.21    
  1.28    
  0.485   
  -0.0377 
  0.419

Specify model parameters

[ param ]  
  WT   = 70
  EGFR = 90
  ALB  = 4.5
  AGE  = 35

[ theta ] @annotated
  0.443   : log(TVKA)
  4.12    : log(TVV2)
  1.17    : log(TVCL)
  4.21    : log(TVV3)
  1.28    : log(TVQ)
  0.485   : log(CLEGFR)
  -0.0377 : log(CLAGE)
  0.419   : log(CLALB)
[ param ]  
  WT   = 70
  EGFR = 90
  ALB  = 4.5
  AGE  = 35
  THETA1 = 0.443   
  THETA2 = 4.12    
  THETA3 = 1.17    
  THETA4 = 4.21    
  THETA5 = 1.28    
  THETA6 = 0.485   
  THETA7 = -0.0377 
  THETA8 = 0.419

Update model parameters

  • Prior to the simulation.
    • param()
    • update()
  • During the simulation

param()

mod$WT
[1] 70
mod <- param(mod, WT=80)
mod$WT
[1] 80
mrgsim(mod) %>% plot("WT")

update()

mod <- update(mod, param = list(WT = 60))
mod$WT
[1] 60
out <- mrgsim_df(mod, param = list(WT = 65))
unique(out$WT)
[1] 65

Update using a object

p <- list(WT = 70.2, FOO = 1) # A list object

mod <- param(mod, p)

mod$WT
[1] 70.2
data <- data.frame(WT = c(70, 80.1), BAR = 2) # A data frame object

mod <- param(mod, data[2,])

mod$WT
[1] 80.1

Update during the simulation

withr::with_seed(12131, 
                 data <- data.frame(
                   ID=1:3, TIME=0, AMT=100, CMT=1, EVID=1, 
                   WT = rnorm(3, 70, 10), EGFR = rnorm(3, 90, 10),
                   ALB = rnorm(3, 4.5, 2), AGE = 50)
)
data
  ID TIME AMT CMT EVID       WT      EGFR        ALB AGE
1  1    0 100   1    1 85.22736  95.47835 1.72194263  50
2  2    0 100   1    1 42.17078 102.73307 4.24730145  50
3  3    0 100   1    1 73.57190  95.16679 0.06696794  50

data_set()

out <- 
  mod %>% 
  data_set(data) %>%
  zero_re() %>% 
  mrgsim(delta = 0.1, end = 24)

plot(out, "WT,EGFR,ALB,AGE,IPRED")

Check if the names match

inventory(mod, data)
Warning: Missing parameters in 'data'
 - THETA1
 - THETA2
 - THETA3
 - THETA4
 - THETA5
 - THETA6
 - THETA7
 - THETA8

[ input ]

New in mrgsolve 1.2.0, you can use the $INPUT block.

[ param ]  
  WT   = 70
  EGFR = 90
  ALB  = 4.5
  AGE  = 35
  THETA1 = 0.443   
  THETA2 = 4.12    
  THETA3 = 1.17    
  THETA4 = 4.21    
  THETA5 = 1.28    
  THETA6 = 0.485   
  THETA7 = -0.0377 
  THETA8 = 0.419   
[ input ]  
  WT   = 70
  EGFR = 90
  ALB  = 4.5
  AGE  = 35

[ param ]
  THETA1 = 0.443   
  THETA2 = 4.12    
  THETA3 = 1.17    
  THETA4 = 4.21    
  THETA5 = 1.28    
  THETA6 = 0.485   
  THETA7 = -0.0377 
  THETA8 = 0.419

Input tags

check_data_names()

{r} # #| echo: true # #| warning: true # check_data_names(data, mod) # check_data_names(data, mod2) # data2 <- rename(data, eGFR = EGFR) # check_data_names(data2, mod2) #

Covariate model


[ prob ]
  ECP8506 mrgsolve model for illustration
  This model requires mrgsolve >= 1.0.3

[ pkmodel ] cmt = "GUT,CENT,PERIPH", depot = TRUE

[ param ]  
  WT   = 70
  EGFR = 90
  ALB  = 4.5
  AGE  = 35

[ theta ] @annotated
  0.443   : log(TVKA)
  4.12    : log(TVV2)
  1.17    : log(TVCL)
  4.21    : log(TVV3)
  1.28    : log(TVQ)
  0.485   : log(CLEGFR)
  -0.0377 : log(CLAGE)
  0.419   : log(CLALB)

[ omega ] @block 
  0.219
  0.0668 0.0824
  0.121  0.0703 0.114

[ sigma ] @block
  0.0399

[ pk ] 
  double V2WT   = log(WT/70.0);
  double CLWT   = log(WT/70.0)*0.75;
  double CLEGFR = log(EGFR/90.0)*THETA(6);
  double CLAGE  = log(AGE/35.0)*THETA(7);
  double V3WT   = log(WT/70.0);
  double QWT    = log(WT/70.0)*0.75;
  double CLALB  = log(ALB/4.5)*THETA(8);

  double KA  = exp(THETA(1) + ETA(1));
  double V2  = exp(THETA(2) + V2WT + ETA(2));
  double CL  = exp(THETA(3) + CLWT + CLEGFR + CLAGE + CLALB + ETA(3));
  double V3  = exp(THETA(4) + V3WT);
  double Q   = exp(THETA(5) + QWT);

  double S2 = V2/1000.0; //; dose in mcg, conc in mcg/mL

[ error ] 
  double IPRED = CENT/S2;
  double Y = IPRED * (1+EPS(1));

[ capture ] CL V2 IPRED Y WT EGFR ALB AGE

C++ functions

if(a == 2) b = 2;
if(b <= 2) {
  c=3;
} else {
  c=4;
}
d = a==2 ? 50 : 100;

C++ functions

double d = pow(base,exponent);
double e = exp(3);
double f = fabs(-4);
double g = sqrt(5);
double h = log(6);
double i = log10(7);
double j = floor(4.2);
double k = ceil(4.2);

\[ d = base^{exponent} \\ e = e^3 \\ f = |-4| \\ g = \sqrt{5} \\ h = \ln{6} \\ i = \log_{10} 7 \\ j = \lfloor 4.2 \rfloor \\ k = \lceil 4.2 \rceil \]

General rules

$PK

double CL = TVCL * pow(WT/70, 0.75) * exp(ETA(1));
bool cure = false;
  • Any valid c++ code is allowed.
  • Declare variable type (e.g., double for numeric variable; bool for logical variable).
  • Use ; behind each line of syntax.
  • Case sensitive, use exp and pow, not EXP and POW.
  • Lots of help on the web (https://en.cppreference.com/w/cpp).

Caution

Be careful of integer division

double result = 1/2; # 0
double result = 1.0/2.0; # 0.5

Random effects


[ prob ]
  ECP8506 mrgsolve model for illustration
  This model requires mrgsolve >= 1.0.3

[ pkmodel ] cmt = "GUT,CENT,PERIPH", depot = TRUE

[ param ]  
  WT   = 70
  EGFR = 90
  ALB  = 4.5
  AGE  = 35

[ theta ] @annotated
  0.443   : log(TVKA)
  4.12    : log(TVV2)
  1.17    : log(TVCL)
  4.21    : log(TVV3)
  1.28    : log(TVQ)
  0.485   : log(CLEGFR)
  -0.0377 : log(CLAGE)
  0.419   : log(CLALB)

[ omega ] @block 
  0.219
  0.0668 0.0824
  0.121  0.0703 0.114

[ sigma ] @block
  0.0399

[ pk ] 
  double V2WT   = log(WT/70.0);
  double CLWT   = log(WT/70.0)*0.75;
  double CLEGFR = log(EGFR/90.0)*THETA(6);
  double CLAGE  = log(AGE/35.0)*THETA(7);
  double V3WT   = log(WT/70.0);
  double QWT    = log(WT/70.0)*0.75;
  double CLALB  = log(ALB/4.5)*THETA(8);

  double KA  = exp(THETA(1) + ETA(1));
  double V2  = exp(THETA(2) + V2WT + ETA(2));
  double CL  = exp(THETA(3) + CLWT + CLEGFR + CLAGE + CLALB + ETA(3));
  double V3  = exp(THETA(4) + V3WT);
  double Q   = exp(THETA(5) + QWT);

  double S2 = V2/1000.0; //; dose in mcg, conc in mcg/mL

[ error ] 
  double IPRED = CENT/S2;
  double Y = IPRED * (1+EPS(1));

[ capture ] CL V2 IPRED Y WT EGFR ALB AGE

Block matrix

[ omega ] @block 
  0.219
  0.0668 0.0824
  0.121  0.0703 0.114

or

[ omega ] @block 
  0.219 0.0668 0.0824 0.121 0.0703 0.114
  • Diagonal elements are variance.
  • Off-diagonal elements are covariance.

\[ \begin{bmatrix} 0.219 & 0.0668 & 0.121 \\ 0.0668 & 0.0824 & 0.0703 \\ 0.121 & 0.0703 & 0.114 \\ \end{bmatrix} \]

Diagonal matrix

[ omega ] 
  0.219 0.0824 0.114

\[ \begin{bmatrix} 0.219 & 0 & 0 \\ 0 & 0.0824 & 0 \\ 0 & 0 & 0.114 \\ \end{bmatrix} \]

@correlation

[ omega ] @correlation
  0.219
  0.0668 0.0824
  0.121  0.0703 0.114
  • Diagonal elements are variance.
  • Off-diagonal elements are correlation, not covariance.

[ omega ] and [ sigma ]

...
[ omega ] @block 
  0.219
  0.0668 0.0824
  0.121  0.0703 0.114

[ sigma ] @block
  0.0399
...
[ pk ]
...
  double KA  = exp(THETA(1) + ETA(1));
  double V2  = exp(THETA(2) + V2WT + ETA(2));
  double CL  = exp(THETA(3) + CLWT + CLEGFR + CLAGE + CLALB + ETA(3));
...
[ error ] 
  double IPRED = CENT/S2;
  double Y = IPRED * (1+EPS(1));

@label

...
[ omega ] @block @labels ETA_KA ETA_V2 ETA_CL 
  0.219
  0.0668 0.0824
  0.121  0.0703 0.114

[ sigma ] @block
  0.0399
...
[ pk ]
...
  double KA  = exp(THETA(1) + ETA_KA);
  double V2  = exp(THETA(2) + V2WT + ETA_V2);
  double CL  = exp(THETA(3) + CLWT + CLEGFR + CLAGE + CLALB + ETA_CL);
...
[ error ] 
  double IPRED = CENT/S2;
  double Y = IPRED * (1+EPS(1));

@annotated

...
[ omega ] @annotated @block
  0.219 : ETA on KA
  0.0668 0.0824 : ETA on V2
  0.121  0.0703 0.114 : ETA on CL

[ sigma ] @block
  0.0399
...
[ pk ]
...
  double KA  = exp(THETA(1) + + ETA(1));
  double V2  = exp(THETA(2) + V2WT + + ETA(2));
  double CL  = exp(THETA(3) + CLWT + CLEGFR + CLAGE + CLALB + + ETA(3));
...
[ error ] 
  double IPRED = CENT/S2;
  double Y = IPRED * (1+EPS(1));

Separated [ omega ]/[ sigma ]

...
[ omega ] 
  0.219

[ omega ]
  0.0824
  
[ omega ]
  0.114

[ sigma ] @block
  0.0399
...
[ pk ]
...
  double KA  = exp(THETA(1) + + ETA(1));
  double V2  = exp(THETA(2) + V2WT + + ETA(2));
  double CL  = exp(THETA(3) + CLWT + CLEGFR + CLAGE + CLALB + + ETA(3));
...
[ error ] 
  double IPRED = CENT/S2;
  double Y = IPRED * (1+EPS(1));

@name

...
[ omega ] @ name KA
  0.219

[ omega ] @ name V2
  0.0824
  
[ omega ] @ name CL
  0.114

[ sigma ] @block
  0.0399
...
[ pk ]
...
  double KA  = exp(THETA(1) + + ETA(1));
  double V2  = exp(THETA(2) + V2WT + + ETA(2));
  double CL  = exp(THETA(3) + CLWT + CLEGFR + CLAGE + CLALB + + ETA(3));
...
[ error ] 
  double IPRED = CENT/S2;
  double Y = IPRED * (1+EPS(1));

Import NONMEM estimates?

  • Use [ nmext ] or [ nmxml ]
    • [ nmext ] reads from the .ext file.
      • Can be faster.
      • Doesn’t retain [ omega ] and [ sigma ] structure.
    • [ nmxml ] reads from the .xml file.
      • Can be slower.
      • Does retain [ omega ] and [ sigma ] structure.

[ nmext ] and [ nmxml ]

If you have a NONMEM script with two $OMEGA

...
$OMEGA BLOCK(2)
0.1 
0.01 0.1

$OMEGA BLOCK(3)
0.1 
0.01 0.1
0.01 0.01 0.1
...
  • [ nmext ] loads OMEGA matrix as a 5x5 matrix.
  • [ nmxml ] loads OMEGA matrix as a 2x2 matrix and a 3x3 matrix.

[ nmext ]

[ theta ] @annotated
  0.443   : log(TVKA)
  4.12    : log(TVV2)
  1.17    : log(TVCL)
  4.21    : log(TVV3)
  1.28    : log(TVQ)
  0.485   : log(CLEGFR)
  -0.0377 : log(CLAGE)
  0.419   : log(CLALB)

[ omega ] @block 
  0.219
  0.0668 0.0824
  0.121  0.0703 0.114

[ sigma ] @block
  0.0399
[ nmext ]  
path = "../nm-model/106/106.ext"
root = "cppfile"

root = "cppfile" look for .ext file relative to where the model source code file is located

[ nmext ]  
run = 106
project = "../nm-model/"
root = "cppfile"

This syntax will let mrgsolve model look for ../nm-model/106/106.ext

[ nmxml ]

[ theta ] @annotated
  0.443   : log(TVKA)
  4.12    : log(TVV2)
  1.17    : log(TVCL)
  4.21    : log(TVV3)
  1.28    : log(TVQ)
  0.485   : log(CLEGFR)
  -0.0377 : log(CLAGE)
  0.419   : log(CLALB)

[ omega ] @block 
  0.219
  0.0668 0.0824
  0.121  0.0703 0.114

[ sigma ] @block
  0.0399
[ nmxml ]  
path = "../nm-model/106/106.xml"
root = "cppfile"

root = "cppfile" look for .xml file relative to where the model source code file is located

[ nmxml ]  
run = 106
project = "../nm-model/"
root = "cppfile"

This syntax will let mrgsolve model look for ../nm-model/106/106.xml

Models in closed form

Like NONMEM, mrgsolve also has “short-cut” to solve one- and two-compartment models with first order input in closed form using [ pkmodel ], which usually results in substantial speed up.

[ pkmodel ]

mrgsolve

[ pkmodel ] cmt = "CENT", depot = FALSE     // 1-cmt IV
[ pkmodel ] cmt = "CENT", depot = TRUE      // 1-cmt oral/SC
[ pkmodel ] cmt = "GUT,CENT", depot = FALSE // 2-cmt IV
[ pkmodel ] cmt = "GUT,CENT", depot = TRUE  // 2-cmt oral/SC

NONMEM

$SUBROUTINE ADVAN1 TRANS2; 1-CMT IV
$SUBROUTINE ADVAN2 TRANS2; 1-CMT oral/SC
$SUBROUTINE ADVAN3 TRANS4; 2-CMT IV
$SUBROUTINE ADVAN4 TRANS4; 2-CMT oral/SC

[ plugin ]

  • autodec
  • nm-vars
  • Rcpp

autodec

Typical way…

[ pk ] 
  double V2WT   = log(WT/70.0);
  double CLWT   = log(WT/70.0)*0.75;
  double CLEGFR = log(EGFR/90.0)*THETA(6);
  double CLAGE  = log(AGE/35.0)*THETA(7);
  double V3WT   = log(WT/70.0);
  double QWT    = log(WT/70.0)*0.75;
  double CLALB  = log(ALB/4.5)*THETA(8);

  double KA  = exp(THETA(1) + ETA(1));
  double V2  = exp(THETA(2) + V2WT + ETA(2));
  double CL  = exp(THETA(3) + CLWT + CLEGFR + CLAGE + CLALB + ETA(3));
  double V3  = exp(THETA(4) + V3WT);
  double Q   = exp(THETA(5) + QWT);

  double S2 = V2/1000.0; //; dose in mcg, conc in mcg/mL

With autodec

[ plugin ] autodec

...

[ pk ] 
  V2WT   = log(WT/70.0);
  CLWT   = log(WT/70.0)*0.75;
  CLEGFR = log(EGFR/90.0)*THETA(6);
  CLAGE  = log(AGE/35.0)*THETA(7);
  V3WT   = log(WT/70.0);
  QWT    = log(WT/70.0)*0.75;
  CLALB  = log(ALB/4.5)*THETA(8);

  KA  = exp(THETA(1) + ETA(1));
  V2  = exp(THETA(2) + V2WT + ETA(2));
  CL  = exp(THETA(3) + CLWT + CLEGFR + CLAGE + CLALB + ETA(3));
  V3  = exp(THETA(4) + V3WT);
  Q   = exp(THETA(5) + QWT);

  S2 = V2/1000.0; //; dose in mcg, conc in mcg/mL

autodec

When this plugin is invoked, mrgsolve will search your model code for assignments and automatically declare them as double precision numbers (double). The following blocks are searched

  • $PREAMBLE
  • $MAIN (or $PK)
  • $ODE (or $DES)
  • $TABLE (or $ERROR)
  • $PRED

Skip autodec a variable

  • When you are using the autodec plugin, you can still declare variables as double or int or bool. mrgsolve already finds those variables and will understand to leave those declarations alone.
[ plugin ] autodec

...

[ pk ] 
  V2  = exp(THETA(2) + V2WT + ETA(2));
  int n = 2; 
  bool ddi = TRUE;
  • In case mrgsolve does try to declare (as double) a variable that shouldn’t be handled that way, you can note this name in an environment variable inside your model called MRGSOLVE_AUTODEC_SKIP.
[ env ] MRGSOLVE_AUTODEC_SKIP = c("my_variable_1")

nm-vars

NONMEM

$PK
...
KA   = EXP(THETA(1)+ETA(1))
V2   = EXP(THETA(2)+V2WT+ETA(2))
CL   = EXP(THETA(3)+CLWT+CLEGFR+CLAGE+CLALB+ETA(3))
V3   = EXP(THETA(4)+V3WT)
Q    = EXP(THETA(5)+QWT) 
...
F1    = LOG(THETA(10))
ALAG1 = SQRT(THETA(11))
...
$DES
DADT(1) = - KA * A(1)
DADT(2) = KA * A(1) - (CL/V)*A(2) - Q*(A(2)/V2 - A(3)/V3)
DADT(3) = Q*(A(2)/V2 - A(3)/V3)
...

mrgsolve

[ pk ]
...
double KA  = exp(THETA(1) + ETA(1));
double V2  = exp(THETA(2) + V2WT + ETA(2));
double CL  = exp(THETA(3) + CLWT + CLEGFR + CLAGE + CLALB + ETA(3));
double V3  = exp(THETA(4) + V3WT);
double Q   = exp(THETA(5) + QWT);
...
F_CENT    = log(THETA(10));
ALAG_CENT = sqrt(THETA(11));
...
[ ode ]
dxdt_DEPOT = - KA * DEPOT; 
dxdt_CENT = KA * DEPOT - (CL/V)*CENT - Q*(CENT/V2 - PERIPH/V3);
dxdt_PERIPH = Q*(CENT/V2 - PERIPH/V3);
...

nm-vars

NONMEM

$PK
...
KA   = EXP(THETA(1)+ETA(1))
V2   = EXP(THETA(2)+V2WT+ETA(2))
CL   = EXP(THETA(3)+CLWT+CLEGFR+CLAGE+CLALB+ETA(3))
V3   = EXP(THETA(4)+V3WT)
Q    = EXP(THETA(5)+QWT) 
...
F1    = LOG(THETA(10))
ALAG1 = SQRT(THETA(11))
...
$DES
DADT(1) = - KA * A(1)
DADT(2) = KA * A(1) - (CL/V)*A(2) - Q*(A(2)/V2 - A(3)/V3)
DADT(3) = Q*(A(2)/V2 - A(3)/V3)
...

mrgsolve with nm-vars

[ plugin ] nm-vars

...

[ pk ]
...
double KA  = EXP(THETA(1) + ETA(1));
double V2  = EXP(THETA(2) + V2WT + ETA(2));
double CL  = EXP(THETA(3) + CLWT + CLEGFR + CLAGE + CLALB + ETA(3));
double V3  = EXP(THETA(4) + V3WT);
double Q   = EXP(THETA(5) + QWT);
...
F1    = LOG(THETA(10));
ALAG1 = SQRT(THETA(11));
...
[ ode ]
DADT(1) = - KA * A(1);
DADT(2) = KA * A(1) - (CL/V)*A(2) - Q*(A(2)/V2 - A(3)/V3);
DADT(3) = Q*(A(2)/V2 - A(3)/V3);
...

nm-vars

  • F, R, D, ALAG
    • To set bioavailability for the nth compartment, use Fn.
    • To set the infusion rate for the nth compartment, use Rn.
    • To set the infusion duration for the nth compartment, use Dn.
    • To set the lag time for the nth compartment, use ALAGn.
  • A, A_0, DADT.
    • To refer to the amount in the nth compartment, use A(n).
    • To refer to the initial amount in the nth compartment, use A_0(n).
    • To refer to the differential equation for the nth compartment, use DADT(n).
  • Math.
    • EXP(a) gets mapped to exp(a).
    • LOG(a) gets mapped to log(a).
    • SQRT(a) gets mapped to sqrt(a).
  • Use T in $DES to refer to the current time in the ode solver rather than SOLVERTIME

Rcpp

  • This gives you functions and data structures that you’re used to using in R, but they work in c++.

  • The main use for this is random number generation. Any d/q/p/r function in R will be available; arguments are the same, but omit n (you always get just one draw when calling from c++.

  • Note: this will slightly increase compile time.

$PLUGIN Rcpp

$ERROR
double u = R::runif(0, 1);

Many other [ plugin ]

  • tad
  • evtools
  • CXX11
  • mrgx
  • RcppArmadillo
  • BH

Resource: https://mrgsolve.org/user-guide/plugins.html

Other blocks

  • [ setup ]
  • [ env ]
  • [ pred ]
  • [ premable ]
  • [ global ]

[ set ]

  • To configure the model object on load.
[ set ] end = 240, delta = 0.5, req = "IPRED"

[ env ]

  • To define a set of R objects that might be evaluated in other model blocks.
  • This block is all R code (just as you would code in a stand-alone R script.
...
$ENV

Sigma <- cmat(1,0.6,2)

rescale <- 1/1000 

convert <- function(x) x * 5 / 166.2
...
$PARAM
p = 500 * rescale
q  = convert(2.51)
...
$OMEGA @object Sigma
...

[ pred ]

  • Use $PRED to write a model without differential equations. In this block, write all algebraic expressions for derived parameters, the response, and any other derived output quantities.

  • It is an error to include the following blocks when $PRED is being used: $MAIN, $TABLE, $PKMODEL, $ODE, $CMT, $INIT.

  • An example NONMEM model without differential equation


$PROBLEM A NONMEM MODEL WITHOUT DES
$INPUT   ID TIME DV AUC
$DATA    SOMEDATA
$PRED
 EVE0=THETA(1)
 E0=EVE0*EXP(ETA(1))
 IMAX=40
 AUC50=100
 RESP=E0-IMAX*AUC/(AUC50+AUC)
 Y=RESP+ERR(1)
$THETA   100
$OMEGA   0.1
$SIGMA   4
$ESTIMATION
  • An example mrgsolve model without differential equation
[ PARAM ] TVE0 = 100, AUC50 = 100, IMAX = 40, AUC = 0

[ PRED ]
double E0 = EVE0*exp(ETA(1));

double RESP = E0 - IMAX*AUC/(AUC50+AUC);

[ premable ]

  • Only called one time during the simulation run (right at the start). The code in this block is typically used to configure or initialize C++ variables or data structures that were declared in [ global ].

[ PLUGIN ] Rcpp

[ GLOBAL ]
namespace{
  Rcpp::NumericVector x;
}

[ PREAMBLE ]
x.push_back(1);
x.push_back(2);
x.push_back(3);

[ MAIN ]
<some code that uses x vector>

[ global ]

  • To define variables outside of any other block
  • Two more-common uses:
    • Write #define preprocessor statements.
    • Define global variables (e.g., double, bool, int)
[ global ] #define CP (CENT/S2)
[ global ] double CP = 0; 
...
[ error ] 
    CP = CENT/S2;

Validation

  • Once we finished developing a popPK model using NONMEM, it is common that we translate a NONMEM model into a mrgsolve model for simulation.

  • Validation needs to be done to confirm the correctness of the translation by comparing:

    • PRED: NONMEM vs mrgsolve.
    • IPRED: NONMEM vs mrgsolve.